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ABSTRACT 


This  report  presents  a statistical  version  of  the  Gaussian 
initial  orbit  technique.  It  neglects  neither  the  angular 
velocity  nor  the  radial  velocity  terras  in  the  f and  g series. 
More  importantly  it  provides  a rigorous,  analytically  simple 
result  for  the  radius  of  convergence  of  the  f and  g series. 

The  radius  of  convergence  can  be  extremely  small  at  periastron, 
approaching  zero  as  the  eccentricity  of  the  orbit  approaches 
unity.  The  leading  term  is  given  by  /2(1  - e)3/|,2P/(3tr)  as 
e -*■  1 where  P is  the  orbital  period.  This  implies  that 
initial  orbit  determination  by  any  procedure  which  uses  the 
f and  g series  is  a process  fraught  with  the  possibility  of 
unknown  errors.  The  central  result  is  given  in  Eq.  (29). 
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I. 


INTRODUCTION 


This  report's  primary  contribution  is  a rigorous  result 
concerning  the  radius  of  convergence  of  the  f and  g series 
of  celestial  mechanics.  It  is  given  in  Eq.  (29).  The  result 
appears  to  be  unknown  in  that  a thorough  perusal  of  the 
classical  texts  (Brouwer  and  Clemence,  Brown,  Brown  and  Shook, 
Danby,  Escobal,  Moulton,  Plummer,  Smart,  Wintner,  etc.)  and 
several  other  texts  reveals  no  mention  of  it.  The  simplicity 
of  the  formula,  its  apparent  connection  with  the  well-known 
divergence  at  e = 0.6627434  associated  with  the  use  of 
Lagrange's  formula  on  Kepler's  equation,  and  its  importance 
lead  me  to  believe  that  1)  there  must  be  a simpler  method  of 
deriving  it  and  2)  while  the  work  here  is  original,  it  must 
duplicate  someone  else's. 

The  f and  g series  are  important  in  orbit  determination, 
the  differential  correction  of  orbits,  and  in  the  numerical 
integration  of  orbits.  Their  use  in  orbit  determination  is 
principally  to  obtain  a distance  estimate  from  the  measurement 
of  angles.  In  the  other  two  cases  the  f and  g series  form 
the  basis  for  the  algorithms  used  in  the  numerical  integration 
of  the  equations  of  motion. 

The  main  thrust  of  the  analytical  result  is  that  for 
highly  eccentric  orbits  (e  > 0.5)  the  radii  of  convergence 
of  the  f and  g series  become  very  small  near  the  periastron 
point.  For  example,  with  an  eccentricity  of  1//2  the  common 


radius  of  convergence  is  [£n(l  + /2)  - 1//2]P/(2tt)  ^ P/36 
(P  is  the  orbital  period) . This  is  ^20m  for  a 2 rev/day 
artificial  satellite.  Hence,  the  usuable  range  of  the  series, 
that  is  when  they  provide  1%  accuracy,  is  <5m  (see  Table  3) . 
Thus,  if  we  remember  to  include  the  inevitable  errors  of 
observation,  the  procedure  of  Gaussian  initial  orbit  determina- 
tion will  be  a difficult  one.  The  neglect  of  the  non-negligible 
velocity  terms  only  compounds  the  problem*. 

As  we  rarely  ever  get  exactly  three  measurements  of  angles 
only,  a statistical  initial  orbit  determination  procedure,  based 
on  the  logic  of  Gauss,  seems  appropriate.  This  is  developed  in 
Sections  II  - V.  It  is  no  longer  necessary  to  throw  away  any 
velocity  effects  using  this  formulation.  The  motivation  for  and 
proof  of  the  central  result  is  given  in  Section  VI.  A brief 
discussion  of  its  implications  (SVII)  and  alternatives  to  the 
Gaussian  technique  (SVIII)  are  also  presented. 

*For  a Molniya  type  satellite  Targument  of  perigee  = 270°, 
inclination  = 60° , eccentricity  = 1//2)  observed  at  the 
equator,  a.  = n/2 , 6 = n/6 , and  t/r  = 2n  where  n is  the  mean 
motion.  Clearly  the  angular  speed  and  the  foreshortening  term 
are  comparable  and  neither  is  small. 
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II.  THE  MOTION  OCCURS  IN  A PLANE 

Let  r(t)  be  the  geocentric  location  vector  in  the  usual 
(approximately)  inertial  reference  frame  where  the  z axis  is 
the  axis  of  the  earth  (e.g.,  the  North  Celestial  Pole  is  at 
x “ 0,  y = 0,  z = + >*>)  , the  z = 0 plane  is  the  extension  of 
the  earth's  equatorial  plane,  the  positive  x axis  points  in 
the  direction  of  the  Vernal  Equinox,  and  the  y axis  completes 
a right  handed  triple; 

r = (x,y,z)  = r(cos5cosa,  cos6sina,  sin6),  (1) 

where  a and  6 are  geocentric  right  ascension  and  declination. 

If  G is  the  constant  of  gravitation  and  is  the  mass  of  the 
earth,  then  the  equations  of  motion  are 

d2r/dt2  = r = - GM0r/|r|3  e - pr/r3  (2) 

where 

r2  = x2  + y2  + z2  = | r | 2 , r > 0.  (3) 
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For  r > 0 there  exists  a unique,  continuous  solution  of  Eqs.  (2) 

for  which  r(t)  takes  on  the  value  r(tQ)  and  v(t)  = r(t)  = dr(t)/dt 

takes  on  the  value  v(t  ).  Moreover,  the  solution  is  a con- 

— o 

tinuous  function  of  £(t0)  and  v(tQ)  and  a continous  function  of 
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small  changes  in  the  right-hand  side  of  Eq.  (2) . These  state- 
ments are  easily  proved  by  using  the  standard  existence  and  unique- 
ness theorems  (§3. Iff  of  reference  1)  for  ordinary  differential 
equations  obtained  from  Picard's  method  of  successive  approxima- 
tions. The  continuity  of  the  solution  with  respect  to  changes  in 
the  force  is  the  foundation  of  perturbation  theory.  The  dependence 
of  the  solution  on  six  independent,  initial  conditions  is  the 
foundation  for  the  statement,  "Three  different  observations  of 
angles  (e.g.,  a and  6)  only  suffice  to  uniquely  determine  an 
orbit" . 

If  we  take  the  vector  cross  product  of  Eq.  (2)  with  r(t) 
we  f ind 

r (t)  X f (t)  = 0,  (4) 


but. 


d[r (t)  X v(t)]/dt  = r (t)  X r (t) , (5) 

so 

L = r ( t)  X v ( t ) , L = | L | = [pa(l  - e2)]1/2,  (6) 

is  a constant  vector.  If  neither  r(t)  nor  v(t)  is  null*,  their 

vector  cross  product  determines  a plane  (reference  3,  831). 

*The  vector  r(t)  can  vanish  only  if  the  motion  is  rectilinear. 
In  that  case  r(t)  and  v(t)  are  collinear,  L vanishes,  and  it's 
obviously  impossible  to  determine  a unique  plane.  If  r > 0, 
then  v(t)  is  never  null.  See  also  reference  2,  §241ff. 
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L is  normal  to  this  plane  and  both  r(t)  and  v(t)  always  lie  in 
this  plane.  Hence,  the  motion  takes  place  in  this  pi which 
is  called  the  invariable  plane.  We  recognize  L as  the  angular 
momentum  (per  unit  mass) . 

Now  consider  three  location  vectors  at  the  three  different 

times  t . , t . , and  t..  Let  r = r(t  ).  Then,  since  (reference  3, 
i j a n — n 

S36)  the  necessary  and  sufficient  condition  that  three  non-null 
vectors,  a,  b,  and  c,  be  coplanar  is 

a * (b  X c)  = 0,  (7) 


it  follows  that 


Ei*  (£j  X r^)  = 0,  i ^ j,  i / k,  j / k.  (8a) 


It  also  follows  from  the  properties  of  the  triple  scalar 
product  (§37  of  reference  3)  that  any  permutation  of  the 
indices  does  not  affect  the  result.  Another  form  of  Eg. (8)  is 


xi 

*i 

z . 

l 

X . 

3 

yj 

z . 
3 

Xk 

*k 
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Since  this  is  a third  order  determinant  there  exists  6 ( = 3!) 
different  ways  of  expanding  it.  In  order  to  demonstrate 
three  of  these  mathematically  redundant  ways  of  expressing 
(9) , define 


Pij  (a,b)  = a Lb.  - ajbi,  (b,a)  = - P.. ^ (a,b)  = Pj:L(a,b),  (10) 


c...  = three  dimensional  Levi-Civita  symbol*, 

1 1 K 


then,  since  D = 0,  three  of  these  ways  are 


° ■ * eijk  *i  pjk(','z> 


(12a) 


= s cijk  yi  !>jklz’x) 


(12b) 


= £ cijk  zi  pjk(x,yJ 


where  i,  j,  and  k are  each  summed  over  1,  2,  and  3.  Even  more 
compact ly , 


(12c) 


D = ? £EabcEijkaiP3k<b'c> 


e123  “ e231  “ G 312  “ 1#  e 132  ' e213  = e321  “lf  a11  other 

= 0 for  i,  j,  and  k one  of  1,  2,  or  3. 


where  E , is  the  alphabetical  counterpart  to  e. ..  and  a,  b,  and 
dbc  13^ 

c are  each  summed  over  x,  y,  and  2. 


In  order  to  remove  the  mathematical  redundancy  of  Eqs.  (12) 
we  employ  the  physics  of  the  problem  to  evaluate  the  P's  and  the 
observations  of  the  problem  to  evaluate  the  remaining  geocentric 
coordinate.  To  do  this  we  need  to  be  able  to  express  r(tn)  in 
terms  of  r(tQ)  and  v(tQ)  at  some  arbitrary  epoch  tQ.  For  this 
purpose  the  f and  g series  are  used. 

Before  proceeding  to  a definition  of  the  f and  g series,  we 

I 

should  note  that  there  is  another  formulation  of  the  coplanarity 
of  the  three  location  vectors  which  is  useful.  Since  the  vectors 
are  three  dimensional,  their  coplanarity  implies  their  linear  | 

dependence.  Thus,  there  must  exist  scalar  constants  d. , d.,  and  f 

1 3 i 

d^  such  that 

P 

(8b) 

« 

!’ 

This  formula  serves  as  the  basis  for  the  Gibbsian  variation  of 
the  Gaussian  technique. 


di^i 


+ d . r . + d.  r.  = 0 . 
D-D  k-k 


III.  THE  f AND  g SERIES 


We  already  know  that,  as  long  as  r > 0,  r(t)  is  a continuous 
function  of  r (t^)  and  v(tQ) . We  choose  to  express  this  in  the 
form 


r(t)  = f[r(tQ),  v(tQ),  t - tQ]  r(tQ) 

+ g[r(tQ) , v(tQ) , t - tQ]v(to) . (14) 

If  we  can  determine  f and  g then  we  can  compute  the  P’s.  The 
usual  method  of  doing  this  is  to  substitute  power  series 
expansions  for  f and  g in  T = t - t into  Eq.  (14)  and  then  use 
that  result  in  Eq.  (2) . Alternatively,  we  note  that 

a = faQ  + gaQ  a = x,  y,  or  z (15a) 

a = - h(r)a  h = y/r3  (15b) 

and  that  all  of  the  higher  derivatives  can  be  computed  by 
successive  differentiation  of  the  equations  of  motion. 

Through  sixth  order  in  T the  result  is 


f - 1 - hQT2/2  ♦ hoPoT3/3  + ho(3qo  - 15p2  ♦ ho>T4/24 

+ hoPo(?Po  " 3<Jo  - ho)T5/8  - V94^  - 210hoPo 

- 630p2qo  + 4 5q2  + 24hoqo  + h2)T6/720,  (16a) 

9 “ T - hQT3/6  + hopQT4/4  + hQ(9qo  + hQ  - 45p2)T5/120 

+ hopo(14p2  - 6qo  - ho)T6/24,  (16b) 

where. 


P - £ * v/r2,  (17a) 

q * v * v/r2  - h.  (17b) 

The  important  thing  to  notice  is  that  the  actual  dependence 
of  f and  g on  r and  v is  really  on  |r  I,  Iv  I , and  |r  * v I. 
The  important  question  to  ask  is  "Do  these  series  converge,  and 
if  they  do,  under  what  circumstances?"  We'll  defer  an  answer 
to  this  question  so  that  we  can  continue  the  discussion  of 
initial  orbit  estimation  (cf.  SVI). 

From  Eq.  (15a)  it  follows  that 
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P..(a,b)  = (f.gj  - Vj)  (aobQ  - aobQ)  = P. . (f ,g)  <aQbo  - aobo) 


" Pij(f^)EabcLc* 


(18) 


Hence,  Eq.  (13)  can  be  written 


D ‘ H Eabccijkaipjk(f-9>Ea  - °‘ 
This  puts  the  physics  into  the  P's.  Since 

r = R + R, 


(19) 


(20) 


where  R is  the  topocentric  location  of  the  celestial  object  and 
R is  the  geocentric  location  of  the  observer,  we  can  replace  a 
by  A + A and  Eq.  (19)  becomes 

D - H Eabc£ijkPjk(t-9>La(Ai  + V ‘ »•  <21> 

Of  the  three  quantities  which  make  up  R^ , two  are  measured.  If 

no  component  of  L vanishes  (i  ^ 0,  ir/2;  ft  / 0,n/2,  it,  3n/2)  , since 

the  P's  are  determined  by  the  physics,  Eqs.  (21)  appear  to  be 

three  linear,  homogeneous  equations  in  the  three  unknowns,  R^ , 

R.,  and  R.  . This,  though,  isn't  true.  We  need  to  examine  the 
3 * 

topocentric  geometry  more  thoroughly  to  explicitly  see  this. 
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IV.  THE  TOPOCENTRIC  GEOMETRY 

At  any  time  t = t^  we  can  write 

rk  = Rk  + r2  + 2RkRcosZk' 
cosZk  = sini{>'sinAk  + cos<J>  'cosAkcos  (ik  - Ak)  , 


(22a) 

(22b) 


where 


R = R (cosAcosA,cosAsinA,  sinA) , (23) 

R = R (cosi)  'cost  ,cos()  ' sinx , sin<)> ' ) . (24) 

Here  xk  is  the  sidereal  time  corresponding  to  tk,  4>*  is  the 
observer's  geocentric  latitude,  and  R is  the  observer's  geocentric 
distance.  There  is  only  one  positive  root  of  Eqs . (22)  and  it 
is  given  by 

Rk  = - RcosZk  + [rk  - R2sin2ZR] 1/2 . (25) 

Therefore,  we  can  replace  the  object's  topocentric  distance  at 
any  time  in  terms  of  known  quantities  and  its  geocentric  distance 
at  that  time.  We  replace  rk  by  its  f and  g series  so  Eqs.  (21) 
explicitly  depend  only  on  rQ,  vq  = 1^1*  and  VQ  = 1 ' Vq|  . 
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The  actual  application  of  Gauss's  technique  involves  the 
explicit  solution  of  Eqs.  (21)  for  one  of  the  R^,  the  use  of 
Eq.  (25) , and  the  total  neglect  of  the  velocity  terms  in  the 
f and  g series.  Thus,  one  equation  in  three  unknowns  is  reduced 
to  one  equation  in  one  unknown.  If  we  actually  look  at  the 
equations  involved,  one  can  (for  nearly  circular  orbits  only) 
heuristically  argue  the  radial  velocity  term  away.  In  any  case 
it's  clear  that  throwing  away  all  of  the  velocity  information  is 
improper.  It's  also  clear  that  the  Gaussian  technique  doesn't 
"solve"  the  initial  orbit  problem. 

Gibbs's  method  starts  from  Eq.  (8b).  He  includes  terms  of 
the  fourth  order  in  time  in  the  f and  g series.  From  Eqs.  (16) 
we  see  that  this  is  also  the  minimum  order  which  includes  all 
the  effects  of  the  velocity.  Thus,  the  Gibbs's  technique  is 
midway  in  complexity  and  accuracy  between  Gauss's  method  and 
the  statistical  method  proposed  in  the  next  section. 
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V. 


STATISTICAL  ORBIT  DETERMINATION 


As  a practical  matter  one  almost  never  acquires  exactly  three 
measurements  of  angles  only.  Therefore,  instead  of  throwing  away 
the  velocity  information,  we  can  solve  for  r =|r|,v  =|v|, 

O 1 “O  1 o 1 — o 1 

and  VQ  = 1^  ' r^|.  To  be  more  explicit  let  us  write  out  Eqs.  (21) 
in  full  for  three,  different,  arbitrary  times  t - t.,t.,  and  t.  ; 

1 j K 


Xijk  = Pjk(f'g)Xi  + Pki(f'g)Xj  + Pij(f'g)xk  = °' 


Yijk  s pjk(f'g)yi  + pki(f>g)yj  + pij(f'g)yk  = °- 


zijk  H pjk(f'g)zi  + pki(f>g)zj  + pij(f'g)zk  = °* 


(26a) 


(26b) 


(26c) 


The  explicitly  appearing  components  of  r . , r . , and  r.  are  to  be 

1 ] K 

replaced  by 


x.  = R.cosA.cosA.  + Rcosd) ' cost  . , 
1111  l 


yi  = RicosAisinAi  + Rcos$ 1 sinx^ , 


Zi  = Ris^nAi  + Rsinf,  etc. 


(27a) 


(27b) 


(27c) 


Each  value  of  R^  is  an  implicit  function  of  r^  given  by  Eqs.  (22) 
or  Eqs.  (22b,  25).  Each  r^  is  given  in  terms  of  rQ,  vq,  VQ,  and 


- tQ  through  the  f and  g series.  We  form  the  quantity 


Ti  = 


S, 


S 


N 

- I 

i, j,k=l 


2 2 

+ YT  ..  + ZT  ..  ) 

ijk  ljk 


(28) 


The  i,j,k  sums  are  unrestricted  from  1 to  N (=  the  total 
number  of  observations).  There  are  3 = N (N  - 1) (N  - 2)/2 
non- zero  terms  in  the  sum.  We  now  minimize  S with  respect  to 
r , v . and  VQ.  Having  found  these  values,  we  then  solve 
for  the  orbital  elements  using  all  of  this  information. 

This  procedure  still  rests  on  the  rapid  convergence  of  the  f 
and  g series  but  no  other  approximation  has  been  made. 

An  efficient  method  of  searching  for  the  minimum  value 
of  S is  the  method  of  steepest  descent.  This  requires  that 
all  of  the  second  derivatives  of  f and  g with  respect  to  h, 
p,  and  q be  computed.  Since  f and  g are  polynomials  in 
these  variables,  this  is  a simple  matter  once  f and  g series  are 
known.  However,  from  the  pattern  in  Eqs.  (16),  one  must 
include  terms  of  the  tenth  order  in  T before  this  can  be 
done  with  any  accuracy.  The  programming  and  testing  of 
statistical  orbit  determination,  for  a wide  variety  of 
different  artificial  satellite  orbits,  is  currently  underway. 

All  terms  inclusive  of  those  twelth  order  in  T are  being  used. 


VI.  THE  RADIUS  OF  CONVERGENCE  OF  THE  f AND  g SERIES 

This  section  has  been  partitioned  into  three  different 
subsections.  We  first  show  (5VIA)  why  we  suspect  that  there 
is  a problem  with  the  f and  g series  for  large  eccentricities. 

We  then  sum  the  f and  g series  in  terms  of  the  true  anomaly 
(I VIB) . Next  (SVIC) , we  provide  a rigorous  proof  for  the 
following : 

The  f and  g series  defined  by  Eq.  (14)  have  the  following 

radius  of  convergence  in  T = t - t ; if  e = 0 the  radius  of 

convergence  is  infinite,  if  e is  unity  the  radius  of  convergence 
3 1/2 

is  (8Q  / 9u ) where  \i  = GMffi  and  Q is  the  distance  from  the 

focus  of  the  parabola  to  its  directrix,  if  ee(0,l)  the  radius 
of  convergence  is  given  by  P*/(2tt)  where  P is  the  period  and 


* 


(29) 


Here  Mq  is  the  value  of  the  mean  anomaly  corresponding  to 
t = tQ,  MQe  . 

In  the  next  section  (§VII)  we  discuss  the  implications  of 
this  for  orbit  determination  using  the  Gaussian  method. 

A.  The  Problem 

We  obtained  the  first  few  terms  of  the  f and  g series  in 
Eqs.  (16).  The  coefficients  of  T are  algebraic  combinations  of 


the  auxiliary  variables  h,  p,  and  q.  In  order  to  get  a 
quantitative  feeling  for  the  relative  size  of  these  three 
quantities  let  us  compute  their  average  over  an  entire  orbit*. 
We  find, 


<h>  = n/[a3(l  - e2)3/2] , (30a) 

<p>  = 0,  <p2>  = e2<h>/2 , (30b) 

<q>  = (1  - e2 ) <h> . (30c) 

Hence,  the  coefficient  of  T in  f is  of  the  order,  on  the  average 

(given  a very  heuristic  averaging  procedure)  of  hk//2.  For  g, 

k ^-1 

the  coefficient  of  T is  of  the  order  of  h 2 . A simple  proof 
by  induction  coupled  with  the  exact  relationship 

dh/dt  = -3hp  (31) 


shows  that  it's  true  for  all  values  of  k.  Hence,  for  the  series 
to  converge,  we  would  need  (by  Cauchy's  root  rest,  §17-4  of 
reference  4) 

T(<h>)1/2<1  or  T/P<  (1-e2)  3/4/(2tt)  . (32) 

Clearly  as  e-*-l  T/P-*0 . 

*For  any  quantity  u,  <u>  = Jrudt/P. 

0 
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The  above  can  be  easily  criticized  especially  for  the  use 

4 

of  averaging  and  the  limit  as  the  eccentricity  approaches  unity. 

The  proper  limit  of  the  formulas  of  elliptical  motion  involves 
both  letting  e-*-l  and  a-*-“>  such  that  Q = a(l-e)  is  finite.  Hence, 
an  asymptotic  expansion  in  1-e  that  starts  from  the  parabolic 
formulas  is  much  more  appropriate.  This  will  be  investigated 
after  we've  summed  the  f and  g series. 

More  generally  though,  despite  our  averaging,  it  appears 
that  the  convergence  of  the  f and  g series  depends  (since  h = p/r3) , 
on  the  satellite's  distance,  e.g.,  Eq.  (32)  without  averaging 
implies  that 


r > 


(T2U>1/3 


(33) 


for  convergence.  This  leads  us  to  ask  the  question,  "What 
fraction  of  an  orbit  is  the  distance  greater  than  some  pre- 
selected lower  bound?".  Suppose  we  parametrize  the  lower  bound 
by  ap  with  p ^ 0.  Then  the  resulting  fraction  is  a function 
of  p and  e only,  say  H(p,e).  A straightforward  computation  yields 


H(P,0) 


1 P < 1, 
0 p > 1. 


H(P,1) 


0 p < Q, 
.1  p > Q. 


(34a) 


(34b) 


j 


1 


I 


” 


\ 


H(p,e)  for  0<e<l  = 


1 


- e . 


- (1/tt  ) { cos-1  [(l-p)/e]  - [e2 


- (1  - p )2r  \ l-e<p<l+e. 


(34c) 


p>l  + e. 


We  note  H(l,e)  = (2e  + tt)/(2tt)  and  3H/3p  < 0 for  ee(0,l) 


H(p,e)  is  tabulated  in  Table  1 for  e = 0(0. 1)0. 9,  p = l-e(0.1)l+e. 


B.  Summing  the  f and  g Series 


The  parametric  representation  of  the  geocentric  position 


x/r  = cosScosa  = cosficosu  - sinftcosisinu. 


(35a) 


y/r  = cos6sina  = sinftcosu  + cosflcosisinu. 


(35b) 


z/r  = sin6  = sinisinu. 


(35c) 


r = a(l  - e )/(l  + ecosv) , 


( 35d) 


u = v + 0), 


(35e) 


tan (v/2)  = [(1  + e) / (1  - e)] 1/2tan (E/2) , 


(35f ) 


E - esinE  = M = n ( t - t ) +M, 

o o 


(35g) 


TABLE  1 

The  Percentage  of  An  Orbit  the  Distance  Exceeds  p • (semi-major  axis) 


The  standard  origin  of  time  is  the  time  of  perigee  passage 
(usually  denoted  by  T) . The  remaining  notation  is  a = semi- 
major axis,  e = eccentricity,  i = inclination,  = longitude 
of  the  ascending  node,  oo  = argument  of  perigee,  v = true 

anomaly,  E = eccentric  anomaly,  M = mean  anomaly,  n = mean 

2 3 

motion,  nP  = 2tt  where  P is  the  period,  and  n a = w. 

If  we  express  Eq.  (14)  explicitly  in  terms  of  the 
parametric  representation  of  the  motion  then  we  find 
(T  = t - t here) 


f = j - r(T)fosin[v(T)  - vj  + r (T)  r^cos  [v  (T)  - vjJ/L,  (36a) 
| r (T)rQsin£v(T)  - vQjj/L,  (36b) 


where  L was  defined  in  Eq.  (6).  All  v's  refer  now  to  the  true 

anomaly,  not  |dr/dt|  and  the  subscript  o indicates  a quantity 

evaluated  at  t = t (e.g.,  v = 0 if  t = T) . Equations 

o o o 

(36)  are  exact  as  long  as  the  Jacobian  3 (x0»y0»z0'*0»^0»^0)/ 

9(a,e,i,ft,uj,Mo)  doesn't  vanish.  Although  we  haven't  computed 

this  determinant,  we  expect  it  to  vanish  for  e = 0,  i = 0, 

i = tt/2.  If  the  motion  is  assumed  to  be  circular  then 

3 2 

p = q = 0,  h = u/a  = n and  by  induction  on  the  coefficients 
of  Tn  we  find 

f = cos[v(T)  - vq] 


# 


(37a) 


r 


g = (1/n) sin [v (T)  - vq] . 


But  if  e = 0,  v=E=M  and  since  the  sine  and  cosine  are 
analytic  for  all  real  values  of  their  argument,  the  power 


series  expansions  of  f and  g in  terms  of  t - tQ  converge 


everywhere.  This  proves  the  first  part  of  the  theorem. 

Let  us  now  consider  the  case  e = 1.  The  analog  of 
Kepler's  equation  is 


1/2 


where  v(tQ)  = 0.  The  explicit  solution  for  tan(v/2)  is 


tan (v/2) 


2 2 ^21 

cV)  J 


= [cT  + (1  + 

2 2 ^/21 
+ cT  - (1  + c TZ)  J 


1/3 


1/3 


with 


c2  = 9y/8Q3,c  > 0,  T - t - tQ. 


(37b) 


tan (v/2)  + (1/3) tan3 (v/2)  = (y/2Q3)  (t  - tQ)  (38) 


(39a) 


(39b) 


. 
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From  this  result,  we  can  compute  sinv  and  cosv  in  terms  of 
T.  If  we  substitute  this  into  Eqs . (36)  and  then  expand  the 
binomials,  we  will  have  the  appropriate  power  series.  They  will 
converge  when  the  binomial  series  converges,  e.g. , for 


c | T | £ 1.  Now  we  can  address  the  earlier  criticism  concerning 
the  limit  as  e-*l.  Since 

Q = lim  a(l  - e) , Q finite,  (40) 

a-*-°° 
e+l 

it  will  be  asymptotically  correct  to  replace  Q by  all  - e)  in  c. 
Thus,  as  e-+l  we  can  expect  convergence  of  the  f and  g series 
when 

3tt  | T | /P  < /2(1  - e)  3/2  • (41) 

Although  of  a different  functional  form  than  the  result  in 
(32),  the  problem  with  e+1  does  not  go  away.* 

Having  provided  this  much  motivation  for  the  existence  of 
a problem,  we  now  turn  to  solving  it  rigorously.  It  is  simpler 
to  work  with  the  eccentric  anaomaly  than  the  mean  anomaly  so 
we  rewrite  the  f and  g sums  as 


f = (a/LrQ)  (rof-Qsinvo  + LcosvQ)  (cosE  - e) 


f 2 1/2  1 , \ 

+ l_a(l  - e ) /r0LJ  vhsinvQ  - rofocosvQ/ sinE, 


(42a) 


*This  is  the  limit  of  Eq.  (29)  when  M = 0 and  e -*■  l" 

o 
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g 


= - (ar  /L)sinv  (cosE  - e) 
o o 


+ ar 

L o 


2 1/2  I 

(1  - e ) /Lj  cosv  sinE. 


(421 


We  need  not  only  to  expand  the  trigonometric  functions  of  the 
eccentric  anomaly  in  a power  series  (which  converges  for  all 
real  values  of  E) , but  to  expand  E in  a convergent  power  series 
in  M,  substitute  this,  term  by  term,  in  the  series  for  cosE  and 
sinE,  rearrange  terms  to  obtain  a power  series  for  f and  g in 
M,  and  then  show  convergence. 

C.  The  Series  f (M) , g (M)  and  Their  Convergence 

The  first  thing  we  need  to  do  is  obtain  E (M)  in  a power 
series.  From  Kepler's  equation,  Eq.  (35g) , we  can  solve  for 
E (M)  in  a convergent  Fourier  series  (§17.2  of  reference  5), 


E 


00 

M + l (2/v) J (ve)sin(vM) , 
v=l  v 


(43) 


where  J^(z)  is  the  Bessel  function  of  the  first  kind,  of  order 
p,  argument  z.  For  e real,  this  series  converges  on  [0,1).  We 

want,  however,  a Taylor's  series  representation  for  E(M).  To 

l l 

compute  the  derivatives  d E/dM  we  differentiate  (43)  l times 

l i . 

term  by  term  and  sum.  The  sum  will  be  d E/dM  if  the  series 
obtained  by  (£  - l)-fold  differentiation  term  by  term  converges 
and  if  the  series  obtained  by  £-fold  differentiation  term  by 
term  converges  uniformly  (§18.5  of  reference  4).  Thus,  we  need 


l 


I 


P 


to  show 


l v*'_1Jv  (ve)  §&§  (vM)  (44) 

v=l 

converges  uniformly  VS.  1 and  ee(0,l)  . The  e = 0,  e = 1 cases 
we’ve  already  disposed  of. 

By  Weierstrass ’ s comparison  test  (reference  6,  §3.34),  since 
for  all  real  M,  |sinvM|  <_  1,  |cosvM|  <_  1,  if 

°° 

£ |v4  J (ve)  | S,  >_  1,  ec(0,l),  (45) 

v=l  v 

converges  then  the  series  (44)  converges  uniformly.  But  Jv(ve) 
for  ee(0,l]  is  a positive  decreasing  function  of  v (reference  5 
§8.5).  It,  therefore,  follows  from  D’Alembert’s  ratio  test 
(reference  4 §17.4)  that  the  series  (45)  converges  whence  the 
series  (44)  all  converge  uniformly  for  S,  1,  ec(0,l].  Therefore 

S.-1  ^ 

+ 2 (-1)  2 l (ve)  cos  (vM)  , S.  odd,  (46a) 

v=l  v 

S, 

— oo 

l)2  l v*'"1Jv(ve)sin(vM)  , S,  even.  (46b) 

v=l 

Since  all  the  derivatives  of  E with  respect  to  M exist,  E 
is  an  analytic  function  of  M (reference  7 §169H  of  Vol.  II).  It 
is,  therefore,  expressible  in  a Taylor's  series  and  the  Taylor 
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series,  in  a sufficiently  small  neighborhood,  does  converge  to  E. 
At  M = 0*  we  have. 


d£E/dMfc | 


M=0  1 


£-1 


6#1  + 2 (-1)  *■  l v£_1J  (ve),  l odd, 
*,A  v=l  v 


0 , l even . 


(47a) 


(47b) 


The  first  few  values  are** 


E ' (0)  = 

1/ (1-e)  , 

(48a) 

E'  ’ ’ (0) 

= -e/ (1-e) 4 , 

(48b) 

EV(0)  = 

e (1  + 9e) / (1-e) 7 , 

(48c) 

EVli(0) 

= -e (1  + 54e  + 225e2) /(1-e) 10. 

(48d) 

Thus , 


E = M/(l-e)  + l Sov  ,M2k  + 1/T(2k  + 2),  (49) 

k=l  1 


where  T(z)  is  the  gamma  function  and 

*The  expansion  about  M = 0 is  equivalent  to  expanding  about  the 
instant  of  perigee  pasSage.  For  any  other  value  of  M we  would 
find  that  the  f and  g series  in  (M  - M ) has  radius  or  convergence 
K.  Hence,  once  the  M = 0 result  is  established  the  MQ  ? 0 
result  follows  immediately. 

**The  reader  can  see  in  Eqs.  (48)  the  kernel  of  the  problem. 
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(50a) 


S2k  + 1 = 2(_1)kv^1Y2k(v,e) 
Y2k(v'e)  E y2kjv (ve) * 


(50b) 


The  radius  of  convergence  of  the  E(M)  power  series,  n,  is 
given  by  (reference  4 §19.5). 

1 A = lim  sup | S_.  ,/r (2k  + 2)|1/(2k  + 1}.  (51) 

k-*»  + x 

Since  for  ee(0,l]  (ve)  is  a positive  decreasing  function  of 

2k 

v and  for  k > 0,  v ^ 1 v is  a positive  increasing  function 
of  v,  it  follows  that  Y2k(v/e)  has  a s^n9le  maximum  as  a 
function  of  v.  Call  this  value  of  v,  N (not  necessarily  an 
integer) . Let  [N]  be  the  greatest  integer  less  than  or 
equal  to  N.  Then, 

m 

Vhi  < IN],  Y2k(m  - i/e)  < / Y2k(v,e)dv  < y2k(m,e),  (52) 

m- 1 

So,  since  JQ(0)  = 1, 


IN]  — 1 ,Ml  IN] 

l Y2k(m/e)  < JY2k(v/e)dv  < I Y2k(ra'e)*  (53) 

m=0  0 m=0 


Also, 
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IN] +m+l 

Vta  > 1,  Y2k(lN]  + m'e)  > / Y2k(v,e)dv  > Y2k  ( tN] +ro+l  ,e)  , 

IN1+m  (54) 


SO 


I Yovlm.e)  > / YoV(v,e)dv  > £ y,.  (m,e) 

m=  [N]  +1  ■**  IN]  +1  ZK  m=2+  [N] 


(55) 


Combining  the  inequalities  (53,  55)  yields 


® [N] +1 

/ Y2k(v,e)dv  - j ^ y2k(v,e)dv  + Y2kdN],e)  + Y2kdN]  + l,e) 


> I Y2k(m,e) . 


m=0 


(56) 


However , 


N 


Y2k(N*e)  > / Y2k(v*e)dv  > Y2k([N],e), 

In] 


(57a) 


and 


IN]  +1 

Y2k(N,e)  > ^ Y,v(v,e)dv  > y,v([N]  + l,e) , 


N 


2k 


2k 


(57b) 


so 
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(57c) 


[N]+l 

{N]  Y2k(v,e)dv  > y2k(lN],e)  + Y2k([N)+l,e) 


Therefore, 


w CO 

J2k  = (jY2k(v'e)dv  >JQ  Y2k(m'e)  “ I s2k  + Xl / 2' 


(58) 


We  now  have,  for  all  k > 0,  ee(0,l]  an  upper  bound  for 

lS2k  + 1l*  The  next  step  is  the  evaluation  of  the  integral  in 
(58)  . 

We  can  write  (reference  5 §8.5) 


J (ve)  = (1A)  / exp [-vF  (0  , e)  ] d6  , 
v 0 


(59) 


where 


F(6,x)  = £n[0  + (02  - x2sin20)1/2]  - £n(xsin6) 

- cot0(02  - x2sin20)  1//2,  (60) 

and  we  can  prove , 

F(0,x)  > F (0 ,x)  > F (0 , 1)  = 0.  (61) 

Replace  Jy(ve)  in  the  integral  of  I2k  with  this  expression, 
interchange  the  order  of  integration  (this  requires  the  uniform 
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convergence  of  the  improper  integral  I2  [reference  6,  14.44]; 
we  shall  explicitly  show  this  below)  and  make  the  change  of 
variable  m = vF(0,x)  in  the  v integral.  Then, 

71  oo 

2I2k  = Id0/F2k+1  <e  p2kexp(-u)du  > |S2k  -t  1|.  (62) 

The  y integral  is  just  T(2k  + 1)  and  from  the  inequality  on 
F(9,e)  above. 


2T(2k  + 1) 

2k+l  J d6  > 

ttF*k  MO.e)  0 


lS2k  + ll 


The  result  for  x is,  then 


(63) 


1/x  = lim 

k -*■<*> 


2T(2k  + 1) 

5FTT 

r(2k  + 2)FZK  A(0,e) 


1/ (2k+l) 


(64) 


Or, 


x = F(0,e) . (65) 

Let  us  now  look  back  on  what  we've  done.  We've  shown  that 
E(M) , expressed  as  a power  series  in  M about  M * 0,  converges  for 
I M | < F(0,e).  The  one  missing  point  is  the  uniform  convergence 
of  the  I 2^  integral  which  we  now  provide. 
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We  need  to  show  that 


J°°  v2k  ■|(1/tt)  J^expl  - vF  (0  ,e)  ]d0Jdv 


= (1/tt)J  jj°°v2kexp[  -vF  (0 ,e) ] dvld0 . (66) 

0 L 0 J 

This  will  be  true  if  the  inner  integral  on  the  right  hand  side 
of  (66)  converges  uniformly  in  0 for  0e[O,ir].  From  (61)  we  can 
majorize  the  integral  by  the  gamma  function  again  so  that  by  de  la 
Vallee  Poussin's  comparison  test  for  convergence  (§6.12  of 
reference  7,  Vol.  I)  the  inner  integral  on  the  right  hand 
side  of  (66)  does  converge  uniformly. 

The  last  thing  we  need  to  do  is  inquire  into  the  permissi- 
bility of  substituting  one  power  series  in  another,  performing 
the  Cauchy  multiplications,  rearranging,  and  summing.  From  Vol. 
II,  §161  of  reference  7 we  can  see  that  in  the  present  circum- 
stances this  is  permissible  and  the  radius  of  convergence  of 
the  f and  g series  is  precisely  A. . q.e.D. 

One  might  wonder  as  to  the  relationship  between  this  result 
and  the  results,  for  instance  in  §100  of  reference  9,  wherein 
r and  v are  expressed  in  a Fourier  series  of  M.  These  are  known 

to  be  divergent  for  some  values  of  M once  e > 0.6627434.  This 

2 

number  happens  to  be  the  modulus  of  a complex  root  of  F (0,e)  = 1. 
The  connection  arises  because  of  the  nature  of  the  Kapteyn 
series  for  E (M)  and  the  use  of  Lagrange's  formula  (reference  10 


1 


556).  The  series  for  r and  v are  really  power  series  in  e whose 
coefficients  happen  to  be  trigonometric  functions  of  the  mean 
anomaly.  If  we  expanded  sinM  etc.  in  a power  series  in  M,  its 

radius  of  convergence  would  clearly  be  at  most  0.66 . In  fact, 

from  reference  11  S4.3,  it  would  be  much  less.  It  appears  that 
the  theorem  given  in  the  beginning  of  this  Section  is  the  most 

I l 

general  result  one  can  obtain  without  special  arguments  depend- 
ing on  the  values  of  e and  M. 

6 
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VII.  THE  IMPLICATIONS  FOR  ORBITAL  ANALYSIS 


F(0,e)  is  given  in  Table  2 for  e = 0(0. 1)1.0.  We  note  that 
F(0,0.1)/F(0,0.7)  = 11  which  shows  the  dramatic  drop  in  the 
permissible  time  span  as  e + 1.  If  n = 2 rev/day,  e = 0.7,  then 
the  maximum  time  for  which  the  f and  g series  converge  is 
<21  minutes.  We  also  give  in  Table  2 the  maximum  eccentricity  such 
that  the  f and  g series  converge  for  |t|/P  <_  0(0.05)0.5. 

Although  it  is  clear  that  for  Mq  = 0,  K -*■  0 very  rapidly  as 

e -*•  1,  we  shall  rarely  be  so  unfortunate  as  to  observe  a satellite 

at  perigee.  The  two  most  promising  search  patterns  will  find 

satellites  near  6=0.  The  most  common  high  eccentricity 

satellites  can  be  characterized  by  (roughly)  n = 2 rev/day, 

e = 1//2,  w = 3it/2.  Hence,  E = ti/4,  Mq  = (tt-2)/4  and  Pn./(2n) 

= 38?3.  Of  course,  V is  a maximum  now  and  v nearly  a maximum 

o o 

so  that  their  neglect  is  especially  serious. 

It  would  seem  that,  since  we  can't  know  the  orbital  phase  of 
the  initial  observations,  the  best  we  can  do  is  retain  a reason- 
able number  of  terms  in  the  f and  g series,  numerically  investigate 
their  divergence,  and  restrict  the  time  span  of  the  observations 
to  a safe,  small,  duration.  Table  3 contains  the  maximum  time 
duration  we  can  use  and  keep  the  relative  percentage  error  at 
1%  and  0.1%  for  e = 0.7  as  a function  of  Mq  = 0(10)60°  for  the 
twelfth  order  f series.  The  results  for  g are  similar.  When  we 
turn  to  the  rapidity  of  the  convergence  of  the  f and  g series  as 


TABLE  2 


RADIUS  OF  CONVERGENCE  OF  THE  f AND  g SERIES 


* 

e 

ft  = F (0  ,e) 

|T|/P 

e 

* 

0.0 

oo 

o 

1" 

0.1 

2.00 

0.05 

0.589 

0.2 

1.31 

0.10 

0.410 

0.3 

0.920 

0.15 

0.293 

0.4 

0.650 

0.20 

0.212 

0.5 

0.451 

0.25 

0.154 

i 

0.6 

0.299 

0.30 

0.112 

0.7 

0.181 

0135 

0.082 

0.8 

0.0931 

0.40 

0.060 

0.9 

0.0313 

0.45 

0.044 

1.0 

0 

0.50 

0.032 
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a function  of  e,  T,  and  Mq  we  discover  that  the  series 
converge  or  they  don't.  What  I mean  by  this  tautology  is  that 
the  partial  sums  of  the  f and  g series,  for  fixed  e,  T,  and 
Mq,  are  either  constant  as  a function  of  their  order  (and 
essentially  equal  to  f and  g)  or  they  are  meaningless.  Numerical 
values  for  the  third  through  twelfth  order  partial  sums  were 
computed  before  reaching  this  conclusion. 

The  other  side  of  the  problem  is  the  desire  to  have  the 
observations  span  as  large  an  extent  of  the  true  argument  of 
latitude  [e.g.,  u of  Eq.  (35e)]  as  possible.  For  highly  elliptical 
orbits  this  will  be  the  case  only  near  periastron,  just  when  the 
radius  of  convergence  is  its  smallest. 


VIII.  ALTERNATIVES 


While  the  statistical  method  outlined  in  §V  is  clearly 
better  than  the  traditional  Gaussian  method,  it  too  relies 
on  the  rapid  convergence  of  the  f and  g series.  With  angles 
only  data  the  other  classical  methods  are  the  Gibbsian 
variation  of  Gauss's  method  and  the  procedures  of  Laplace 
and  Lagrange.  I briefly  mentioned  the  Gibbsian  variation  of 
Gauss's  method  on  page  12.  Since  it  does  minimally  include  the 
velocity,  it  is  better  than  Gauss’s  procedure.  However,  as  it 
too  relies  on  the  rapid  convergence  of  the  f and  g series,  it 
is  subject  to  the  above  criticism  of  the  statistical  Gaussian 
method.  The  Laplace  method  requires  two  numerical  differentiations 
of  the  topocentric  direction  cosines.  I don't  think  this  is 
a good  thing  to  do.  The  Laplace  method  does  not  rely  on  the 
f and  g series  though.  The  Lagrangian  method  only  uses  a 
single  numerical  differentiation  of  the  topocentric  direction 
cosines  but  does  use  the  f and  g series.  Hence,  this  is 
open  to  the  same  difficulties  as  the  Gaussian,  Gauss-Gibbs, 
and  statistical  Gaussian  techniques. 

If  one  has  angles  and  angular  rates,  which  should  be 

the  case  for  artificial  satellites,  one  can  improve  the 

situation.  A modified  Laplace  method  has  been  developed 

wherein  only  one  numerical  differentiation  is  required 

(reference  12) . A modified  Lagrange  method  without  any 

12 

numerical  differentiations  has  been  developed  but  it  still 


relies  on  the  f and  g series.  Finally,  there  is  yet  another 
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method  which  is  exact  (given  that  the  force  is  -yr/r  ) 

and  only  uses  two  sets  of  observations  of  angles  and  angular 

rates.  Unfortunately,  the  exact  technique  is  very  susceptable 

to  errors  in  the  angular  rates.  [To  make  this  statistical 

would  require  f and  g series.] 

If  we  step  back  from  all  of  this  and  remember  the  history 

of  the  first  three  hundred  years  of  dynamical  astronomy,  we 

realize  almost  all  of  the  development  has  been  concerned 

with  perturbations  of  a circular  orbit  with  the  observer  in 

the  orbital  plane.  As  the  practical  problem  is  to  compute 

initial  orbits  for  highly  eccentric  motion,  it  would  appear 

that  a perturbation  of  a parabolic  orbit  is  the  analysis  to 
9 

pursue.  Moulton  has  made  a start  on  this. 
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NOTE  ADDED  IN  PROOF 


After  the  full  review  of  this  technical  note  Donald  Batman 
showed  me  a reference  in  the  second  edition  of  Moulton's  book. 

It  occurs  under  the  heading  "The  Laplacian  Method  of  Determining 
Orbits"  (pg.  202)  and  refers  to  a paper  of  Moulton's  (published 
after  the  first  edition  of  the  text)  in  which  one  can  find  "the 
determination  of  the  exact  realm  of  convergence"  of  the  f and  g 
series.  The  full  reference  is 

F.  R.  Moulton,  Astron.  J.  2^3,  93  (1903). 

In  this  paper,  by  a method  completely  independent  of  the 
technique  used  herein,  Eq.  (29)  is  derived.  He  also  gave 
tables  similar  to  Table  2 but  was  discussing  orbit  determina- 
tion for  asteroids,  not  artificial  satellites. 
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